////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  -----------         Ebrahim Foulaadvand and Somayyeh Belbasi, 02 June 2012       ---------------      //
//                                                                                                        //
// The programme "AirDragFall" evaluates the impact velocity of a freely fallen object released from a    //
// height h above the ground in the presence of air force in the constant gravitational field g.          //
// The Euler and Euler-Richardson algorithms are implemented. Both the linear and quadratic air drag      //
// forces are simulated. In the acceleration function choose the type of air drag force and inactive      //
// the other type. The paramters are explained in the "FreeFall" programme.                                                                                       //                                                                                                                 //
//                                                                                                        //                                                                                                         //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;

double  tau=0.001,g=9.8,m=1.,c1=0.1,c2=0.001,v0=0.,h,delh=2,Y,vhit,ac=0.,yEU,vEU,yER,vER,tmid=0.,ymid=0.,vmid=0.;

int N=20000,t,s,S=200;

double a (double &v); // acceleration function

main()

{

ofstream filevEU("vEU c2=0.001.plt");

ofstream filevER("vER c2=0.001.plt");

ofstream filevAN("vAN c2=0.001.plt");

ofstream filevFree("vFree.plt");

for (int s=0;s<=S;s++){

h=s*delh;

// ----------- initial conditions --------------------

yEU=h;   vEU=v0;

yER=h;   vER=v0;


//------------------------ Euler Method------------------------

Y=h; t=0;

while (Y>0.01){

yEU=yEU + tau*vEU;

vEU=vEU + tau*a(vEU);

Y=yEU;

t+=1;

} //end while

vhit=vEU;

filevEU<<h<<" "<<fabs(vhit)<<endl;

//cout<<"vhit EU= "<<fabs(vhit)<<endl;//getch();


//------------------------ Euler-Richardson Method------------------------

Y=h; t=0;

while (Y>0.01){

tmid=t*tau + 0.5*tau;
ymid=yER + 0.5*tau*vER;
vmid=vER + 0.5*tau*a(vER);

vER=vER + tau*a(vmid);
yER=yER + tau*vmid;

Y=yER;

t+=1;

} //end while

vhit=vER;

filevER<<h<<" "<<fabs(vhit)<<endl;

double u=exp(h*c2/m);

double chinvU=log(u+sqrt(u*u-1));

filevAN<<h<<" "<<sqrt(m*g/c2)*tanh(chinvU)<<endl;

filevFree<<h<<" "<<sqrt(2*g*h)<<endl;


}// end s

cout<<"end."<<endl; getch();

return 0;

}

//----------------------------------------------------------------

double a (double &v) {


//return ac=-g-(c1/m)*v; // linear air drag

return ac=-g-(c2/m)*fabs(v)*v;  // quadratic air drag

}